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\ ABSTRACT 


Of  particular  importance  in  an  interactive  curve  and  surface  design 
system  is  the  interface  to  the  user.  The  mathematical  model  employed  in 
the  system  must  be  sufficiently  flexible  for  interaction  between  designer 
and  machine  to  converge  to  a satisfactory  result.  The  mathematical  theory 
of  Total  Positivity  is  combined  with  the  interactive  techniques  of  Bezier 
and  Riesenfeld  in  developing  new  methods  of  shape  representation  which 
retain  the  valuable  variation-diminishing  and  convex  hull  properties  of 
Bernstein  and  B-spline  approximation,  while  providing  improvements  in 


the  interactive  interface  to  the  user.  Specifically,  extending  the 
Bezier  notion  of  using  a polygon  to  describe  a smooth  curve,  methods  of 
assigning  a weight  to  each  vertex  which  will  control  the  amount  of  local 
fit  to  the  polygon  or  polygonal  net  are  provided.  Thus,  the  designer  can 
cause  “cusps^and  ‘'Tlats^easily  by  manipulating  the  ’’’tension ""‘at  each 
vertex.  Further,  the  generalization  from  curves  to  surfaces  can  be  done 
with  rectilinear  data  or  triangular  data.  Illustrations  are  provided 
from  an  experimental  implementation  of  the  newly  constructed  models  as 
a demonstration  of  their  feasibility  and  utility  in  computer  aided  curve 


and  surface  design. 

V 
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I.  INTRODUCTION 

Computer  Aided  Geometric  Design  is  principally  concerned  with 
the  modelling  of  physical  objects  within  a digital  computer  for 
automated  design  and  manufacture.  Potentially,  the  computer  can  free 
the  designer  from  the  limitations  of  traditional  drafting  techniques 
while  enhancing  and  accelerating  the  production  process.  Working 
interactively  at  a graphics  terminal  or  a numerically-controlled 
drafting  machine,  a designer  can  modify  or  manipulate  an  existing 
model  or  he  can  design  an  object  ab  initio,  using  trial  and  error 
techniques  to  produce  an  acceptable  design.  The  computer  system 
can  then  generate  the  necessary  information  for  numerically  controlled 
manufacture  andproduction  [1]. 

Statement  of  Problem 

The  two  central  problems  this  paper  addresses  are: 

(1)  The  development  of  a unifying  theory  for  the  construction 
and  analysis  of  methods  of  modelling  free-form  curves  and 
surfaces  in  a digital  computer,  and 

(2)  The  application  of  thie  theory  toward  the  construction 
of  new  techniques  for  ab  initio  design. 

These  problems  are  part  of  a study  termed  Computer  Aided  Geometric 

Design  (CAGD). 
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Computer  Based  Modelling  for  Geometric  Design 

The  realization  of  such  a system  for  curve  and  surface  descrip- 
tion imposes  certain  constraints  on  the  computer-based  model.  The 
model  must  accurately  represent  a variety  of  shapes,  be  amenable  to 
analysis  and  manipulation,  and  must  take  into  account  the  capabilities 
and  limitations  of  both  the  designer  and  the  computer.  Although  it 
is  clear  a mathematical  model  is  required,  the  properties  of  shape 
cannot  be  characterized  entirely  by  the  properties  of  analytic  func- 
tions and,  therefore,  new  mathematical  techniques  for  synthesizing, 
storing  and  retrieving  shape  information  are  needed. 

The  shapes  we  are  concerned  with  in  CAGD  have,  in  general,  no 
specific  functional  form  and  what  is  required  is  an  acceptably  close 
* fit  which  maintains  the  character  of  the  object  being  designed.  Since 

all  real  three-dimensional  objects  are  of  finite  size,  the  mathemati- 
cal representation  of  these  objects  should  reflect  this  boundedness  in 
a convenient  and  efficient  manner.  Shape  is  non-analytic  in  nature, 
in  that  it  is  not  generally  possible  to  predict  the  entire  shape  of 
an  object  from  its  shape  in  any  one  localized  area.  Further,  shape 
is  an  orientation  independent  phenomenon,  thus  intrinsically  a proper- 
ty of  the  object  itself  and  not  of  the  space  (coordinate  system)  in 
which  it  is  embedded. 

To  meet  these  shape  criteria  we  restrict  our  attention  to  vector- 
valued bounded  piecewise  analytic  functions  from  a finite  dimensional 
linear  space.  That  is,  we  have  for  curves 
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P(t)  - l <fr  (t)P., 
iel  i 1 


and  for  surfaces 


(1.1) 


P(u,v)  = Z <k-(u,v)P.,  (1.2) 

id  1 1 

3 

where  I is  some  linear  ordered,  finite  set  of  integers,  P.  e R , t and 

2 

(u,v)  are  elements  of  some  bounded  subset  of  R and  R , respectively, 
and  the  {^(t)}  and  {<J>.(u,v)}  are  linearly  independent  sets  of  bounded 
piecewise  analytic  functions. 

The  set  of  points  P^  clearly  constitutes  the  controlling  parame- 
ters with  respect  to  the  basis  {<Jk,  iel}.  It  remains,  then,  to  find 
bases  which  appropriately  model  the  constraints  of  a CAGD  system.  Not 
only  must  the  control  parameters  be  suitable  for  interactive  communica- 
tion between  designer  and  computer,  but  the  basis  functions  must  be 
appropriate  for  digital  representation  and  computation.  At  this  point 
we  make  the  distinction  between  the  fi tting  of  a mathematical  represen- 
tation to  a predesigned  object  and  the  design  of  an  object  ab  initio, 
for  our  choice  of  basis  may  depend  on  which  type  of  design  process  we 
wish  to  model. 

In  fitting  we  are  concerned  with  finding  a mathematical  description 
of  an  existing  shape  or  physical  object,  and  in  ab  initio  design  the 
problem  is  to  create  a mathematical  representation  which  meets  design 
constraints  which  may  be  entirely  subjective,  esthetic  considerations. 
Because  we  know  the  exact  shape  a priori , fitting  can  be  done  by  digi- 
tizing points  off  the  object  and  interpolating  the  data  with  an 
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appropriate  approximation  method.  What  constitutes  a "good"  approxi- 
mation is  ultimately  up  to  the  designer  and  interaction  with  the 
model  may  be  necessary.  Thus  there  are  hard  and  soft  constraints  to 
be  met.  The  hard  constraints  are  the  discrete  data  which  must  be 
interpolated,  while  the  soft  constraints  are  manipulated  by  the  de- 
signer in  determining  an  acceptable  fit.  In  general,  there  are  few 
hard  constraints  in  ab  initio  design  and  the  designer  manipulates 
this  increased  flexibility  to  achieve  his  final  design. 

While  interpolation  techniques  appear  necessary  for  the  fitting 
problem,  this  requirement  is  not  present  in  ab  initio  design.  Of 
fundamental  importance  in  ab  initio  design  is  the  interface  to  the 
i designer.  If  an  interactive  design  session  is  to  be  successful,  the 

designer  must  be  able  to  manipulate  the  control  points  in  an  itera- 
tive process  which  converges  when  a satisfactory  shape  is  achieved. 

Thus  the  shape  of  the  curve  or  surface  must  respond  predictably  to 
the  interactive  manipulation  of  the  control  points.  Experience  with 
interpolation  techniques  indicates  that  unwanted  kinks  and  oscillations 
can  occur  during  the  design  process,  making  it  difficult  to  predict 
the  response  of  the  curve  to  manipulation  of  the  data. 

An  ab  initio  system  for  the  design  of  automobiles  which  avoids 
these  difficulties  is  Bezier's  Systeme  Unisurf  at  Kenault  [1].  In  Systeme 
Unisurf  the  control  points  are  associated  with  the  vertices  of  a poly- 
gon for  curves  and  the  vertices  of  a rectilinear  net  for  surfaces. 

„ The  resulting  curve  or  surface  replicates  the  gross  features  of  the 

polygon  or  rectilinear  net.  More  importantly,  the  design  responds 


predictably  to  the  maniuplation  of  the  control  points. 

In  Bezier's  method  curves  are  piecewise  polynomials  where  the 
curve  is  given  on  each  piece  by 


m 

p(t)  = z <j>.(t)P.  , (l. 

1=0  1 

where 


4-U)  = (?)  t1  (l-t)m_1  , 

for  te(0,l),m  is  the  degree  of  the  polynomial,  and 'fhe  P.  are  the 
appropriate  vertices  for  that  piece.  Figure  1.1  shows  a typical 
Bezier  curve  and  the  associated  Bernstein  basis  functions. 


Figure  1.1  (a)  Cubic  Bezier  curve  and  associated  polygon. 
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Surfaces  are  generated  by  taking  the  tensor  product  of  the  basis  func 
tions  with  respect  to  two  orthogonal  directions,  i.e., 


P(u,v) 


where  (u,v)e(0,l )x(0,l ) , m,n  are  the  degrees  of  the  polynomial  with 
respect  to  the  u and  v directions  and  the  P..  are  the  vertices  of  the 

• vJ 

associated  rectilinear  network.  Figure  1.2  shows  a Bezier  network  and 
the  corresponding  surface. 


7 


Figure  1.2  Bicubic  Bezier  surface  and  associated  recti- 
linear network  of  points. 

i Gordon  and  Riesenfeld  [14]  have  associated  the  remarkable’  repro- 

ducing power  of  the  Bernstein-Bezier  methods  with: 

(1)  the  fact  that  Bernstein  approximation  is  variation- 
diminishing,  and 

(2)  the  fact  that  the  Bernstein  basis  functions  are  posi- 
tive and  sum  identically  to  1,  i.e., 

m 

l Mx)  = 1 

i=0 

and 

^.(x)  > 0 

for  all  x,  such  that 

xe(0,l)  and  i = 0,1,  ....  m.  (1.5) 
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If  we  let  V(f)  denote  the  number  of  sign  changes  of  some  function 
f and  Vfx^.x^.  ....  xm]  denotes  the  number  of  sign  changes  of  the  in- 
dicated sequence,  then  property  (1)  is  equivalent  to 

V(P)  < V[P0,Pr  ....  PJ  (1.6) 

Property  (2)  is  commonly  referred  to  as  the  convex  hull  property, 
since  the  conditions  stated  are  necessary  and  sufficient  for  any 
curve  of  the  form  (1.1)  to  lie  within  the  convex  hull  of  the  coeffi- 
cients P^ 

Note  that,  in  general,  the  piecewise  Bernstein  approximation  is 
only  C°  continuous.  Although  it  is  not  difficult  to  "fix  up"  con- 
tinuity, many  applications  in  the  automobile  and  aircraft  industries 
2 

? require  at  least  C continuity.  Riesenfeld  [5]  has  recently  proposed 

vector-valued  B-spline  approximation  as  a generalization  of  Bernstein 
approximation  which  retains  the  variation-diminishing  and  convex  hull 
properties  while  improving  the  continuity  class  of  the  curves  and  sur- 
faces to  arbitrary  smoothness. 

The  vector-valued  B-spline  approximation  of  degree  M-l  to  the 
associated  polygon  [pi , » •••»  Pn]  for  integral  knot  spacing  is  given 

by 

P(t)  = " *iM(t)  P,  (1.7) 

i=l  1M  1 

3 

where  P^eR  , te(-°°,«>),  n > M and 
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** 


1 

MI 


E (-1)J  (") 
j=0  J 


(s-j-i) 


M-l 

+ 
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(1.8) 


where 


(s  - (J+i))J 


(s-j-i )M  for  s > (i+j) 
0 otherwise 


(1.9) 


As  with  Bezier's  method,  the  surface  equation  is  given  as  the  tensor 
product 


P(u,v) 


^ *jN(u) 


+iM(v>  P 


„)■ 


(u,v)eR2 


M 

n > N,  m _>  M and  <J>.  represents  the  ith  B-spline  basis  function  of 

degree  (M-l),  as  in  (1.8).  An  important  property  of  the  B-spline 

basis  of  degree  (M-l)  is  that  each  basis  function  is  already  a piece- 

M 2 

wise  polynomial  of  continuity  C , thus  the  designer  does  not,  in 
general,  have  to  explicitly  join  segments  to  form  the  compound  curve 
(surface),  although  the  piecewise  nature  of  the  curve  ( surface ) is 
implicit  in  the  basis.  Figures  1.3  and  1.4  show  typical  B-spline 
curves,  while  Figure  1.5  shows  a B-spline  surface  and  the  associated 
rectilinear  network. 
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We  noted  earlier  that  interpolation  methods  seem  ill-suited  for 
ab  initio  design  in  that  it  may  be  difficult  to  control  interactively 
the  occurrence  of  kinks  and  oscillations  in  the  design,  while  it  ap- 
pears the  variation-diminishing  methods  of  Bezier  and  Riesenfeld  per- 
form well.  The  mathematical  statement  of  these  observations  is  cap- 
tured in. the  following  exclusion  theorem  due  to  Schoenberg  [8]. 

Theorem  1.1.  Let  L(f)  be  a linear  transformation  defined  for  all  con- 
tinuous functions  f(x)  on  [0,1],  L(f)  f f for  some  f.  Then  if 
L ( a + bx  + cx  ) = a + bx  + cx  , for  all  a,  b,  c e R (reals),  then  L 
cannot  be  variation-diminishing. 

From  the  point  of  view  of  ab  initio  design  then,  techniques  of  approx- 
imation which  have  a high  degree  of  "reproductive  power"  cannot  be 
variation-diminishing,  and  thus  there  is  no  guarantee  that  you  can 
"control"  oscillations  during  manipulation  of  the  vertices  of  the 
defining  polygonal  network.  Since  most  polynomial  and  polynomial 
spline  techniques  of  interpolation  are  precise  for  more  than  linear 
polynomials,  i.e.,  they  reproduce  polynomials  of  degree  two  and  higher, 
they  cannot  be  variation-diminishing.  Thus  Schoenberg's  theorem  estab- 
lishes a firm  foundation  for  approach  to  CAGD.  If  the  problem  is  to 
fit  hard  constraints,  then  interpolations  methods  would  seem  superior 
to  variation-diminishing  methods,  while  variation-diminishing  methods 
appear  superior  for  ab  initio  design. 

Although  vector-valued  B-spline  approximation  offers  an  attractive 
generalization  of  Bezier's  methods  for  curve  and  surface  design,  both 
techniques  possess  inherent  weaknesses.  The  extension  from  curves  to 
surfaces  requires  a rectilinear  network  of  points.  Thus  a designer  is 


13 


' 


restricted  to  rectilinear  data  as  control  points.  There  is  no  simple 
way  of  adding  and  deleting  points  within  this  restrictive  topology. 

In  general,  the  designer  must  delete  or  add  a whole  row  and  column 
of  points  to  maintain  the  rectilinear  topology  of  the  data.  This  is 
clearly  unacceptable  for  it  would  generate  global  changes  to  the 
surface.  Further,  Bernstein  and  high  degree  (continuity)  B-spline 
approximation  techniques  have  low  "reproductive"  power  in  the  sense 
that  the  curves  are  poor  approximations  to  the  shape  of  the  defining 
polygon  [5].  The  ability  to  control  interacti vely  the  amount  of 
local  fit  to  the  polygon  while  maintaining  a high  continuity  class 
would  be  an  important  feature  of  an  ab  initio  design  system. 

This  paper  is  an  investigation  into  the  development  of 
i mathematical  models  which  retain  the  valuable  variation-diminishing 

and  convex  hull  properties  of  Bernstein  and  B-spline  approximation, 
while  providing  improvements  in  the  interface  to  the  user  for  ab  initio 
design.  A family  of  prototypes  which  includes  the  Bernstein  and 
Bezier  method  is  presented.  Specifically,  methods  of  assigning  a 
weight  to  each  vertex  which  will  control  the  amount  of  local  fit  to 
the  polygon  or  polygonal  net  are  provided.  Thus  the  designer  can 
easily  cause  "pseudo-cusps"  and  "pseudo-flats"  by  manipulating  the 
"tension"  at  each  vertex.  Further,  the  generalization  from  curves  to 
surfaces  can  be  done  with  rectilinear  data  or  triangular  data.  There- 
fore, the  difficulty  involved  with  the  addition  and  deletion  of  points 
in  a rectilinear  network  can  be  avoided.  The  feasibility  and  utility 
^ of  the  newly  constructed  models  for  interactive  design  are  demonstrated 

on  an  experimental  system  for  curve  and  system  design. 


14 

A unifying  theory,  based  on  total  positivity  [3],  is  developed 
here  for  the  construction  of  such  models.  Total  positivity  is  a con- 
cept that  has  played  an  important  role  in  various  mathematical 
sciences.  The  application  of  the  theory  of  total  positivity  to  CAGD 
has  resulted  not  only  in  new  techniques  for  synthesizing  shape  opera- 
tors for. CAGD  but  also  has  provided  general  methods  for  analyzing  the 
Bernstein-Bezier  and  B-spline  methods  previously  discussed.  Section 
II  provides  a general  discussion  of  the  representative  and  manipula- 
tion rules  for  totally  positive  functions  and  their  interrelationships. 
In  Section  III  we  apply  the  theory  of  total  positivity  to  the  analysis 
of  some  existing  methods  of  ab  initio  design  and  then  proceed  to  the 
construction  of  new  shape  operators  for  curve  and  surface  design.  In 
Section  iv  we  discuss  the  utility  of  these  methods,  as  demonstrated  on 
the  curve  and  surface  design  system. 


II.  TOTAL  POSITIVITY — A REVIEW 


Total  Positivity  and  the  Variation-Diminishing  Property 

In  Section  I we  established  the  usefulness  of  variation-diminish- 
ing methods  for  ab  initio  design.  Due  to  the  importance  of  Schoenberg's 
exclusion  theorem  (1.1)  in  this  connection,  we  restate  and  prove  it 
here. 

Theorem  2.1  [8].  Let  T(f)  be  a linear  transformation  defined  for 

all  f(x)  continuous  on  an  interval  [a,b],  where  g(x)  = T(f)  is  itself 
a continuous  function  in  [a,b],  with  the  following  properties: 

1.  T(c  + dx)  = c + dx,  for  all  c,  d e reals. 

2.  V(T(f))  < V(f). 

3.  There  exists  t such  that  T(t)  t f,  i.e.,  T is  not  the  iden- 
tity transformation. 

Then  there  exist  real  numbers  a,  6,  y such  that 

T(a  + gx  + yx^)  t a + Bx  + yx^, 
for  xe[a,b]. 

Proof:  Assume  the  converse,  that  g(x)  = T(f)  does  preserve  quadratic 

polynomials,  i.e.,  T(x  ) = x . Let  a < xQ  < b.  We  will  show  that 

g(xQ)  = f(xQ),  which  implies  T must  be  the  identity  transformation,  in 

contradiction  of  (3)  above. 

«# 

Indeed,  assume,  for  instance,  that 

* 

- 
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g(x0)  > f(x0) 

since 

2 2 
T(f(x)  - a - bx  - cx  ) = g(x)  - a - bx  - cx 

by  the  variation-diminishing  property  we  have 

V(g(x)  - a - bx  - cx2)  < V(f(x)  - a - bx  - cx2).  (2.1) 

That  is,  the  graph  of  g(x)  crosses  the  graph  of  any  quadratic  poly- 
nomial no  more  often  that  the  graph  of  f(x)  does.  But  for  a such  that 

f(xQ)  < a < g(xQ)  , 

p 

the  parabola  y = a + B(x  - Xq)  enjoys,  for  sufficiently  large  posi- 
tive 6,  the  following  properties: 

(i)  It  is  crossed  by  y = g(x)  at  least  twice. 

(ii)  It  is  not  crossed  by  y = f(x). 

2 

(Consider  the  zeros  of  g(x)  - (a  + B(x  - xQ)  and  f(x)  - (a  + 

2 

B(x  - Xq)  .)  These  conclusions  contradict  the  variation-diminishing 
property  (2.1)  and  the  theorem  is  established. 

We  have  already  alluded  to  the  strong  interrelationship  between 
variation-diminishing  methods  and  totally  positive  functions.  This 
relationship  is  established  in  the  following  theorem  due  to  Karlin  [3]. 
Theorem  2.2  [3].  Let  K(x,y)  be  a function  of  two  variables  x e X 
and  y e Y,  where  X and  Y are  linearly  ordered  subsets  of  the  real  line 
and  consider  the  transformation 

g(x)  * T(f)(x)  = ^ K(x,y)  f(y)  dy,  (2.2) 
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where  f and  K are  bounded  and  Riemann  integrable  on  Y.  When  Y con- 
sists of  a discrete  set,  we  interpret  the  integral  as  a sum. 


If  K is  totally  positive  of  order  r,  then 

V(g)  < V(f)  provided  that  V ( f ) < r-1.  (2.1) 

Note  that  equation  (1.1)  is  a vector-valued  form  of  (2.2)  where  Y is 
a discrete  set.  The  proof  of  Theorem  2.2  follows  as  a corollary  to 
the  statement  that,  for  a totally  positive  matrix  A and  any  vector  x 
(of  compatible  length),  Ax  has  no  more  sign  changes  than  does  7.  Be- 
fore we  can  prove  either  of  these  assertions  a general  discussion  of 
the  representation  and  manipulation  rules  for  totally  positive  func- 
tions is  necessary. 

Total  Positivity 

Definition  2.1.  A real  function  K(x,y)  of  two  variables  ranging 
over  linearly  ordered  sets  X and  Y,  respectively  is  said  to  be  totally 
positive  of  order  r (abbr.  TPr)  if  for  all  sequences  x-j  < x^  < ... 

< xm,  yi  < y^  < •••  < ym,  x.eX,  y^.eY,  1 < m < r we  have  the  inequali- 
ties 


K(xi  ,yi ),  ...,  K(xi , yffl) 

K<vV K(v 


> 0. 

(2.4) 


If  strict  inequality  holds  in  2.4,  we  say  K(x,y)  is  strictly  totally 
positive  of  order  r (abbreviated  STPr ) . The  subscript  is  normally 
omitted  when  a function  is  (strictly)  totally  positive  of  all  orders. 
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Note  that  total  positivity  implies  positivity.  As  examples  of  totally 
positive  functions,  we  have: 

(a)  exy,  X x Y e R2  is  STP 

(b)  K^x.i)  = ('J)xi(l-x)m'i  is  STP  of  order  m+1  on  i e I = 

{0,1 , ....  m)  and  x e (0,1 ) . 

(c) .  The  square  n x n matrix  A = [a ( i , j ) = 6..]  is  TP. 

• J 

The  verification  if  (c)  is  trivial , while  (b)  wil  1 be  shown  to  fol  low 
immediately  from  (a).  Note  that  the  function  K(x,i ) in  (b)  represents  the 
Bernstein  basis  of  degree  m-1,  and  the  proof  of  (b)  in  conjunction 
with  Theorem  2.2  would  be  sufficient  to  show  that  Bernstein-Bezier 
methods  are  variation-diminishing.  We  provide  the  proofs  for  (a)  and 
(b)  below. 

Theorem  2.3.  exy  is  STP  on  X x Y = R2. 

Proof:  We  first  show  that 


eVi 

eW 


exly" 

exnYn 


t 0 


for  x-| , X£*  ...  xn  distinct  in  X and  y^ , y^,  ...,  yn  distinct  in  Y. 
This  fact  follows  immediately  if  we  can  show  that  any  exponential 
polynomial  of  the  form 

" a.eyix  (2.5) 

i=l  1 


where 

n 2 
l a.  > 0 
i=l  1 

has  at  most  n-1  zeros. 
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Denoting  the  number  of  distinct  zeros  of  a function  by  Z(f),  we  prove 
(2.5)  as  a lemma. 

Lemma  2.1.  Let  y^,  i=l,  ...,  n,  be  a set  of  distinct  real  numbers. 


Then 


a/1*)  i"-1 


where 


E a.  > 0 . 
i=l  1 

Proof:  (by  mathematical  induction)  Let  n = 1 . Then  clearly  aieyix 
has  no  zeros. 

Now  assume  the  hypothesis  holds  for  n = k-1 , i.e., 


(S!  a'eV) 


< k-2 


for  all  a^eR  such  that 


k-1  ? 

E a/  > 0 
i=l  1 


We  must  show  that 


(Ji  a<eV) 1 w 


for  all  a^cR  such  that 


E a.  > 0. 
i=l  1 


(2.6) 


(2.7) 
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t 

t 


* 


Assume  without  loss  of  generality  a-j  f 0.  Since  e~yix  has  no  zeros, 
we  have 


k 

E 

i=l 


aieyi 


Differentiating  e"yix 


v x 

E a.eJi  with  respect  to  x,  we  get 
i=l  1 


(2.8) 


where 


and 


k 

E 

i=2 


a!eyi 


I 


a((yry,) 


(2.9) 


which,  by  our  induction  assumption, has  at  most  k-2  zeros.  It  follows 
from  Rolle's  theorem  that  a function  can  have  at  most  one  more  dis- 
tinct zero  than  its  derivative,  therefore, 

Z ^ E aieyiX^  < k-1  (2.10) 

Now  we  show  that  for  < x2  < ...  < xn»  and  y^  < y2  < •••  < 
y , y^eY,  the  determinant 

e(X’ M 

\yr  •••»  yn  ' 

cannot  achieve  both  negative  and  positive  values. 

Lemma  2.2.  Given  x^  < x2  < . . . < xn,  x^X  and  y1  < y2  < ...  < yn 
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y^eY,  where  X x Y = R , we  have 


x, , . . . , x 
1 n 


'r  V 


exlyl  . . . eVn 
exnyl  . . . exnyn 


(2.11) 


is  of  one  strict  sign. 

Proof:  This  will  be  proved  by  contradiction.  Fixing  y^  < y^  < . . . < 
yn*  yjeY’  let  xi  < x2  < < xn’  and  Z1  < z2  < *•*  < zn’  xi’zi  in  X 


be  such  that 


V xn 


'r  •••’  yn 


> 0 > e 


z, , ....  z 
1 n 


yr  — ’ yn 


(2.12) 


Since  ex  is  a continuous  real  function,  we  have  for  some  X e(0,l)  the 
determinant 


rr 


D(X)  = e | j = 0 

Xx^  + (1-X)  Zy  ... , Xxn  + ( 1 - X ) zn 


By  (2. 5) this  is  impossible  for  distinct  x^  and  z..  Therefore,  there 


exist  i,j,  i j , such  that 


Xx.  + (1-X)  z.  = Xxj  + (1-X)  z j 


which  implies 


0 = 


1-X  x^Xj 


contradicting  the  strictly  increasing  nature  of  the  x^  and  z^ 
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n 


I 

i 


0 


Now  all  that  remains  to  complete  the  proof  of  Theorem  2.3  is  to 
exhibit  a set  of  increasing  x.  and  increasing  z^,  for  which  the  deter- 
minant in  (2.11)  is  positive.  If  we  let  x^  = i,  i = 0,  1,  ....  n, 
and  y^  = lnz.  where  z.  is  a strictly  increasing  sequence  of  positive 
reals  numbers,  then  (2.11)  specializes  to 


exoyo  exoyn 

1 z'  . . . z" 

0 0 

i 

eWo  eVn 

= n (z.-z.)  > 0,  (2.13) 

i>j  J 

since  the  determinant  on  the  right  hand  side  is  the  well-known 
Vandemonde  determinant  [15]. 

In  order  to  prove  the  total  positivity  of  the  Bernstein  basis, 
we  need  the  following  two  composition  rules  for  totally  positive  func- 
tions. 

Theorem  2.4  [3].  Let  K(x,y)  be  TPr(STPr)  on  X and  Y and  let  <p(x)  and 
4;(y)  be  nonnegative  (positive)  on  X and  Y,  respectively.  Then 

L(x,y)  = 4>(x)  i|>(y)  K(x,y)  (2.14) 

is  TPr(STPr). 

Proof:  The  conclusion  (2.14)  follows  immediately  from  the  fact  that 
the  determinant  is  a multilinear  function  of  its  rows  or  columns. 
Theorem  2.5  [3].  Let  K(x,y)  be  TPr (STPr ) on  X and  Y and  let  <J>(u)  and 
*(v),  ucU,  veV  be  strictly  increasing  functions  mapping  U and  V into 


X and  Y,  respectively.  Then 
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L(u,v)  = K[<J>(u) » tMv)]  (2.15) 

is  TPr(STPr),  ueU,  vcV. 

Proof:  The  result  (2.15)  follows  immediately  from  the  obvious  fact 
that  <p  and  4*  map  increasing  sequences  into  increasing  sequences. 
Theorem  2.6.  The  Bernstein  basis  K^^x)  = ('I')  x\l-x)m’^ 

i-0,1 m,  xe(0,l)  is  TPm+^  on  {0,1,  ....  m)  and  (0,1). 

Proof:  Rewriting  K^i.x)  = (?)  x’O-x)"-1  as 

•Snd.x)  = (?)  (l-x)m  e1  ]‘x  , (2.16) 

the  conclusion  follows  directly  from  Theorems  2.3,  2.4  and  2.5. 

We  will  find  the  following  composition  formula  for  totally  posi- 
tive functions  of  fundamental  importance  in  the  uevelopment  of  shape 
operators  for  CAGD. 

Theorem  2.7  [3].  Let  K(x,z)  and  L(z,y)  xeX,  ycY  and  zeZ  be  Riemann 
integrable  functions  of  z,  where  X,  Y and  Z are  linearly  ordered  sets 
of  the  real  line,  and  define 

M(x,y)  = / K(x,z)  L(z,y)  dz, 

Z 

where  again  we  interpret  the  integral  as  a sum  when  Z is  discrete. 
Then,  if  K(x,z)  is  TPr  on  X x Z and  L(z,y)  is  TP$  on  Z x Y,  then 
M(x,y)  is  TPt  on  X x Y,  where  t = min(r,s). 

Proof:  The  conclusion  follows  easily  from  the  generalized  Cauchy  - 
Binet  formula  for  the  composition  of  determinants  given  by 
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* 


*# 


where  M,  K and  L are  defined  as  above  and  we  assume  Z = [a,b]  for  simpli- 
city. Formula  (2.17)  is  a direct  extension  of  the  Cauchy-Binet  formula 
[17]  for  matrices,  which  we  restate  below. 

Theorem  2.8  [17].  Let  A,  B and  C denote  matrices  of  real  numbers 
or  orders  n x m,  n x k,  and  k x m,  respectively.  If  A = BC,  then 


(2.18) 


where  the  lefthand  term  stands  for  the  minor  of  A involving  rows  i^.i^, 
....ip  and  j1 ,j2,. . . ,jp,  respectively. 

Proof:  The  conclusion  will  follow  if  we  can  establish  the  correspond- 
ing result  for  any  square  matrix  which  is  the  product  of  two  rectangu- 
lar matrices.  So,  suppose  that  a square  matrix  C = ||c^j||'j1  is  the 
product  of  two  rectangular  matrices  A = ||a^||  and  B = II b^ j ||  of 
dimension  m x n and  n x m,  respectively.  That  is, 
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* 


iJ 


n 

E 

a =1 


a . b . 
on  aj 


(i,j  = 1,2,  ...,  m). 


(2.19) 


By  (2.19)  the  determinant  of  C can  be  represented  in  the  form 


C11  •••  clm 


cml  ‘ ' ' cmm 


n 

E aia  1 ' 

<^=1  |al  V 


n 

* ' ^ , ala  ba  m 

a =1  m m 
m 


n n 

E a„  b , . . . E a_  b 
met,  a,  1 


_i  met,  a,  l 'ma  a m 

a,=l  in  a_=l  m m 


ll 


m 


n 

E 


a,  b , ...  a,  b 
la  a 1 la  am 

11  mm 


3 b , ...  a b 

mci|  Oj  1 my  cyn 


/I  2 ...  m \ 

ba/a, 

\aia2...am/  1 < 


a m 
m 


(2.20) 


If  m > n,  then  among  the  numbers  a-| , a2»  ...»  am  there  are  al- 
ways at  least  two  that  are  equal,  so  that  every  summand  on  the  right- 
hand  side  of  (2.20)  is  zero.  Hence  in  this  case  C = 0. 

Now  let  m $ n.  Then  in  the  sum  on  the  right-hand  side  of  (2.20) 
all  those  summands  will  be  zero  in  which  at  least  two  of  the  subscripts 
0| , a^,  ....  am  are  equal.  All  the  remaining  summands  of  (2.20)  can 
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i 


i 


1 


* 


be  split  into  groups  of  m!  terms  each  by  combining  into  one  group 
those  summands  that  differ  from  each  other  only  in  the  order  of  the 
subscripts  ^so  within  each  such  group  the  sub- 
scripts a i , a have  one  and  the  same  set  of  values).  Now 

within  one  such  group  the  sum  of  the  corresponding  terms  is 


where  k,  < k0  < ...  < k is  the  normal  order  of  the  subscripts  a,, 

1 Z m l 

N 

a,,  ...,  % and  efcy  a,,  ...,  c^)  = (-1)  , where  N is  the  number  of 
transpositions  of  the  indices  needed  to  put  the  permutation  ci|,  a2, 

...  « into  normal  order. 

Hence  (2.19)  implies  (2.18), 

Variation-Diminishing  Transformations 

We  are  now  in  a position  to  prove  Theorem  2.2.  As  noted  earlier, 
Theorem  2.2  is  a direct  consequence  of  the  corresponding  theorem  for 
matrices,  which  we  establish  as  Theorem  2.9  below.  The  proof  given 
here  is  essentially  that  given  by  Schoenberg  [16,21].  In  the  follow- 
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ing  definitions  and  theorems,  sequences  and  matrices  will  often  be 
represented  with  the  functional  notation  x(i),  a(i,j),  emphasizing 
the  fact  that  we  are  dealing  with  functions  over  discrete  sets. 
Definition  2.2.  Given  a function  f(t)  defined  in  I,  an  ordered  set 
of  the  real  line.  We  define  the  variation  of  f over  I as 

V(f)  = sup  V[f (t1 ) , f(t2),  ...  f(tm)]  , 

where  the  supremum  is  extended  over  all  sets  t,  < t0  ...<  t (t.rl), 
m arbitrary  but  finite  and  V[x(l),  x(2),  ....  x(m)]  is  the  number  of 
sign  changes  of  the  indicated  sequence,  zero  terms  discarded.  By 
convention  V[0,0,  ....  0]  = -1.  ' 

Definition  2.3.  Given  a sequence  X = [x ( 1 ) , x(2),  ...,  x(n)]  of  real 
numbers,  we  define  N[x(l),  x(2),  ....  x(n)]  as  the  number  of  x(i)'s 
which  are  nonzero. 

Definition  2.4.  A real  m x n matrix  A = [a ( i , j ) ] is  said  to  be 
variation-diminishing  if 
n 

y(i)  = E a(i,j)  x(j),  x(i  = 1,2,  ....  m) 

j=l 

implies  that 

V[y(l),  y(2),  ....  y(m)]  < V[x(l),  ....  x(n)]  . 

Lemma  2.1.  Consider  the  following  matrix  operations: 

(i)  Multiplication  of  a row  or  column  by  a non-negative  constant; 

(ii)  Addition  of  a row  or  column  to  an  adjacent  row  or  column; 

(iii)  Omission  of  a row  or  column. 

These  operations  when  applied  to  a TP  matrix  yield  a TP  matrix. 


28 

Proof:  (iii)  is  trivial,  while  (i)  follows  from  Theorem  2.4.  ( 1 i ) is 
evident  since  any  minor  in  the  new  matrix  would  be  the  sum  of  minors 

in  the  old. 

Lemma  2.2.  Let  X = [x(i),  i = l,2,...,m]  be  a sequence  of  real  num- 
bers. If  for  all 

i0  < ii  < ...  < ip,  P < m, 

we  have 

V[x(i0),  ....  x(ip)]  < p - 1,  (2.21) 

then 

V[x(l),  x(m) ] < p - 1 . 

Proof:  Assume  the  converse,  that  there  exists 

<o  < V<  •••  < Vl 

such  that 

V[x(iQ),  ....  x(ip+1 )]  > p - 1 

then 

V[x(1Q),  ....  x(ip+1 )]  = p,  (2.22) 

since  all  other  cases  immediately  contradict  (2.21).  But  to  preserve 
(2.21)  for  all  subsequences  of  [x( iQ) , ....  x(ip+1)]  the  sequence  must 
alternate  in  sign.  However,  we  than  have  V[x(iQ),  x(1  f^)]Jp+i 

a contradiction  to  (2.22).  From  here  a simple  induction  argument 
yields  the  desired  result. 
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Lemma  2.3.  Let  A be  a TP,  real  m x n matrix  of  rank  p and  define 


y(i ) = £ a(i,j)  x(j),  i = 1,2,  ....  m 

j=l 


where  X = [x(j),  j = 1,2,  ....  n]  is  some  real  sequence.  Let  A' 
be  the  submatrix  of  A consisting  of  the  rows  iQ,  i^,  ....  ip  and 
define 


ot(  K ) = A 


V •••'  Vr  Vr  •••■  V 


for  come  selection  of  columns  j^,j2,  ...»  jp-  Then 


Proof:  Since 


E (-D  ct(k)  y(i.)  = 0 

k=0  K 


y(ij  = E a(i.  j)  x(j)  , 
K j=l  K 


we  have 


? H)k  <*(k)  y(i.) 
k=0  K 


P l n 

c / 1 \ * _ / I.  \ r» 


E (-1)  a(k)  E a(i.,j)  x(j) 
(=0  j=l  K 


" (-Dj  x(j)  ( E (-1)k+j  «(k)  ati.jA  . 
j=l  ' k=0  k / 


(2.23) 


But  for  j not  a member  of  •••»  jp]  the  inner  sum  is  the 
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expansion  of  some  minor  of  A of  order  p+1  and  thus  is  zero.  Similar- 
ly. for  j a member  of  ....  jp],  the  inner  sum  corresponds  to 

the  expansion  of  the  determinant  of  some  matrix  which  possesses  two 
identical  columns  and  again  the  sum  is  zero.  Therefore,  (2.23)  must 
be  identically  zero. 

Theorem  2.9  [ 21  . If  a real  m x n matrix  A = [ a ( i , j ) j 

is  TP,  then  it  is  variation-diminishing. 

Proof:  We  first  assert  that  it  is  sufficient  to  demonstrate  this 

under  the  assumption  that 

V[x(l ),. . . ,x(n)]  = N[x(l ) , . . . ,x(n)]  - 1 . 


% 


For  if  any  x(k)  = 0 we  compress  X and  A by 


x'(j)  = 


j < k 
j > k 


(j  ~ 1.  • • • » n-1 ) 


(2.24) 


a'O'.j) 


a(i.j)  j < k 
a(i,j+l)  •.  j > k 


(j  = 1,  ....  n-1) 


We  then  have  y(i)  defined  by 

n-1 

y(i)  = I a ' (i , j)  x'(j)  . 
j=l 


The  matrix  [a'(i,j)]  is  TP  by  Lemma  2.1  and  if  we  can  show 


V[y(l),  .... 


y(m)]  < V[x ' (1 ) , .... 


x‘ (n-1 )], 
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then  since 

V[x ’ (1 ) , x’(n-l)]  = V[x(l ) , ....  x(n)] 

it  will  follow  that 

V[y(l),  ....  y(m)]  < V[x(l),  ....  x ( n ) ] . 

So  we  can  assume  no  x^  is  zero. 

Next,  we  can  assume  x(k)  and  x(k+l)  are  of  opposite  sign.  If  this  is 
not  true,  then  there  exists  A > 0 such  that  x(k+l)  = Ax(k).  If  we 
compress  X and  A by 


x'(j) 

x(j) 
x(j+l ) 

j < k 

J > k 

(j  = 1,  ....  n-1) 

1 a(i,j) 

j < k 

(2.25) 

a'(i.j) 

| a ( i , k ) + 

Aa(i ,k+l ) 

j=k  (j~l>  ...9 

n-1) 

, a(i , j+1 ) 

j > k 

then  we  have 

n-1 

y(i)  = 2 

j=i 

a ‘ ( i . j ) 

x'(j)  . 

The  matrix  [a'(i,j)]  is  TP  by  Lemma  2.1  and  if  we  can  show  that 
V[y(l),  ....  y(m)]  < V[x'(l),  ....  x'(n-l)], 

then  since 

V[x'(l),  ....  x'(n-l)]  = V[x(l),  ...,  x(n) ] , 


it  will  follow  that 
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V[y(l),  . ..,  y(m)]  < V[x(l),  ....  x(n) ] 

We  have  shown  that  we  may  assume  (2.21).  We  complete  the  proof  then 
with  the  following  lemma. 

Lemma  2.4.  [21].  If  a real  matrix  A = [a(i,j)],  (i=l,  ....  m, 

j-1 .....  n)  is  TP,  then 

V[y(l),  ....  y(m) ] < N[x(l),  ....  x ( n ) ] - 1 

From  above  we  may  assume  no  x^  is  zero.  That  is,  we  may  assume 
N[x(l),  ...,  x(n)]  = n.  Since  r(A)  (rank  of  A)  < n,  we  would  be  done 
if  we  can  show  that 

V[y(l),  ...,  y(m)]  < r(A)  - 1, 

We  prove  this  result,  in  turn,  by  the  following  lemma. 

Lemma  2.5.  C21 1 - If  a real  m x n matrix  A = l a ( i , j ) 1 

is  TP,  then 

V[y(1),  ....  y(m)]  < r(A)  - 1 

Proof:  We  proceed  by  induction  on  r(A).  The  result  is  clearly  true 
if  r(A)  = 0 or  1.  Suppose  that  it  has  been  established  for  r(A)  = 0, 
1,  ....  P-1.  We  must  show  that  it  is  true  for  p. 

We  may  suppose  that  m > p,  since  if  m = p,  then 
V[y(l ).. . . ,y(m)]  < p-1  trivially.  By  Lemma  2.2 

it  is  enough  to  show  that,  if  1 < iQ  c i,  < ...  < ip  < m,  then 
V[y(ig)»  ....  y(ip)]  < P-1.  If  A'  is  the  submatrix  of  A,  consisting 
of  the  rows  ig,  i^,  ...»  ip,  then  A'  is  TP.  There  are  two  cases: 
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r(A')  < p and  r(A')  = p.  If  r(A')  < p,  then 

V[y(i0) y ( ip) 3 < r(A')  - 1 < p - 1 


by  our  induction  assumption. 

So. we  assume  r(A')  = p.  Let  j,,...,jp  be  a selection  of  columns 
of  A'  and  set 


for  0 £ k £ p.  Since  A1  is  TP,  no  two  a(k)'s  are  of  opposite  sign 
. and  because  r(A')  = p,  we  can  so  choose  jj,...,jp  such  that  not  all 

I a(k)'s  are  zero.  We  have  from  Lemma  2.3 

I «(0)  y ( i o ) - a(  1 ) y ( i -j ) = ...  + (-l)p  a(p)  y( ip)  « 0 

(2.26) 

| 

It  is  easily  seen  that  V[y(iQ),  y(i] ) . . ,y(ip)]  = p is  not  compat- 
ible with  (2.26)  so  that  V[y(iQ) , . . . ,y( ip)]  < p -1 . Thus  we  have 
proved  Lemma  2.5,  thereby  proving  Lemna  2.4,  thereby  proving  the 
theorem. 

■ 

Theorem  2.2  is  established  in  the  following  three  corollaries 
* of  Theorem  2.9. 

Corollary  2.9.1 . Let  [y ( 1 ) ,. . . ,y(m)],  [x(l ),. . . ,x(n)]  and  A be 
| * given  as  in  (2.20)  where  now  we  assume 

i 

1 
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! 


i 


i 


i 


V[x(  1 ) . . ,x(n)]  £ r-1 , r£n-l, 
and  A is  TPr.  Then 

V[y(l ) , . . . ,y(m) ] < V[x(l) x(n)]. 

Proof:  Using  the  techniques  given  in  (2.24)  and  (2.25)  we  can 
compress  A and  [x(l ),. . . ,x(n)]  such  that 

N 

y(i ) = E a'(i,j)  x * ( j ) 
j=l 

where  V[x'(l ),. . . ,x' (N)]  = N-l  = V[x(l ) ,. . . ,x(n)]  and  A'  = [a'(i,j)] 

is  a TP  mxr  matrix.  The  conclusion  now  follows  from  Theorem  2.9. 

n 

Corollary  2.9.2.  Let  f(x)  = I 4>(i.x)  a..  If  (<J>(i,x),  i=l,...,n> 

i=l  1 

is  TP  on  x = [a,b]  and  I = {l,2,...,n},  then  (<J>(i,x),  i=l,...,n} 
is  variation-diminishing. 

Proof:  We  must  show 

V(f (x) ) < V[a1,...,an] 

That  is,  for  any  finite  sequence  of  points  x.,  j=l,...,m  we  must 

vJ 

show 

V[f(x1),...f(xm)]  < V[a1 ,an] 

where 

n 

f ( x j ) = E <j>(i  ,x.)  a. 

J i=1  J 1 


I 
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for  j=l,...,m.  But  this  is  exactly  the  result  of  the  above  theorem 
and  we're  done. 

Corollary  2.9.3.  Let  <t>(x,y)  and  g(y)  be  bounded  and  continuous 
functions  on  [a,b]  x [c,d]  and  [c,d],  respectively.'  Let 

d 

f(x)  = / ( x ,y ) g(y)  dy 

c 

Then,  if  <p(x,y)  is  TP,  we  have  V[f(x)]  £ V[g(y)]. 

Proof:  We  must  show  V[f (x1 ) , . . . ,f (xm) ] <_  V[g (y-j ) , . . . ,g(yn)] 

U x-  e[a,b],  y.  e[c,d],  m and  n fixed  but  arbitrary. 

We  approximate  f(x.)  arbitrarily  closely  by  the  Riemann  sum 

J 

n 

f 1 ( x ^ ) = h E 4>(x^  ,yj)  g(yj), 

J 1 

such  that  sgn(f-j(x.))  = sgn(f(x^))  for  i=l,...,m.  The  conclusion 
now  follows  from  Corollary  2.9.1. 

Corol lary  2.9.3  [3].  If  in  Theorem  2.9  we  have,  in  fact  , 

V[y(l),...,y(m)]  = V[x(l),...,x(n)], 

for  A STP,  then  {y ( i ) , i=l,  i=l,...,m)  and  (x(j),  j-l,...,n}  exhibit 
the  same  arrangement  of  signs. 

Proof:  From  the  above  arguments  we  know  it  is  enough  to  assume 

V[x ( 1 ) , . . . ,x(n)]  = n-1.  Under  this  assumption  we  need  only  show 

that  the  first  component  of  [x(l ) ,. . . ,x(n)]  has  the  same  sign  as  the 

first  nonzero  component  of  [y(l),...,y(m)].  Choose  i,  < ^ < • 

n 

such  that  y.  y.  < 0.  Since  y.  = E a.  • x.  (v  = l,2,...,n), 

\ W+l  \ j=l  vJ  J 
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we  have 


Mow  the  minors  of  A are  strictly  positive  and  the  values  (-l)v+^ 

y.  (v  = l,2,...,m)  maintain  the  same  sign.  Thus  sgn(x.)  = sgn(y.  ). 
v 1 ’l 

Theorem  2.2  combined  with  Theorem  2.6  provide  enough  power  to 

conclude  that  scalar-valued  Bernstein  approximation  is  variation- 
diminishing.  However,  as  we  established  in  Section  I,  we  are 
interested  in  vector-valued  approximation  and  transformation  methods, 
thus  we  must  provide  some  extension  of  the  scalar-valued  theory  to 
vector-valued  mathematics.  This  is  provided  in  the  following  defini- 
tion and  theorem  due  to  W.  W.  Meyer  [5]. 
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Definition  2.5  [5].  A vector-valued  approximation  or  transformation 
method  is  variation-diminishing  if  it  is  variation-diminishing  as  a 
scalar-valued  method. 

Theorem  2.10  [5].  If  a vector- val ued  transformation  is  invariant 

under  euclidean  transformation,  then  no  (hyper)  plane  is  pierced  more 
often  by  a vector-valued  transformation  than  by  the  primitive  curve 
itself. 

Proof:  Surely  the  theorem  is  true  for  any  principal  (hyper)  plane 
x(i)  = 0 because  it  is  variation-diminishing  in  the  particular  coordin- 
ate function  x^t).  But  any  (hyper)  plane  can  be  designated  as  a prin- 
cipal one  by  a euclidean  transformation.  Hence  the  result  is  true  for 
all  (hyper)  planes. 

It  is  not  difficult  to  show  that  transformations  of  the  form  (1.1) 
are  invariant  under  euclidean  transformation  if  and  only  if 

W.x)  = 1 . (2.27) 

i 

This  fact  does  not  form  a formidable  obstacle  to  us,  since,  given  a TP 
basis  <f> ( i , x ) , by  Theorem  2.4  we  can  construct  another  TP  basis 

(2.28) 

(2.29) 


<J>'(i»x)  = <j>(i  ,x)  / 2 $ ( j , x ) 

j 

such  that 

Z <J»'(i,x)  = 1 
i 


We  have  established  the  useful  fact  that  the  Bezier  method  is 
variation-diminishing.  In  Section  III  we  continue  with  the  analysis 
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of  other  mathematical  schemes  in  use  or  proposed  for  use  in  ab  initio 
design.  We  then  apply  the  theory  of  total  positivity  to  the  construc- 
tion of  a whole  family  of  new  models  for  CAGD. 


III.  APPLICATION  OF  TOTAL  POSITIVITY  TO  CAGD 

Through  the  application  of  the  theory  of  total  positivity  we  have 
shown  in  Section  n that  Bernstein-Bezier  methods  are  variation-dimin- 
ishing. From  Section  I we  know  that  widely  used  techniques,  such  as 
polynomial  and  spline  interpolation,  are  not  variation-diminishing.  In 
this  chapter  we  continue  this  type  of  analysis  by  applying  the  tech- 
niques developed  in  Section  II  to  other  methods  in  use,  or  proposed 
for  use,  in  aj)  initio  design. 

Piecewise  Bernstein-Bezier  Methods 

In  general,  it  is  vector-valued  piecewise  Bernstein  approximation, 
not  Bernstein  approximation,  which  is  used  in  CAGD.  As  noted  in 
Section  i,  shape  is  essentially  non-analytic  in  nature,  therefore,  it 
is  not  surprising  that  we  often  find,  when  working  interactively,  that 
a curve  segment  is  not  sufficiently  flexible  to  adopt  a desired  shape- 
We  may  then  either  increase  the  order  of  approximati or.  (and  thus  that 

of  the  polygon)  or  the  segment  may  be  split  into  two  or  more  segments, 
(see  Figure  3. 1 ).  Curve  splitting  is  simple  mathematically  and  has  the 
advantage  computationally  of  retaining  a reasonable  (polynomial)  order 
of  approximation.  Since  the  resulting  curve  is  a piecewise  Bernstein 
approximation,  the  question  arises  whether  this  basis  is  variation- 
diminishing  as  well.  We  establish  this  fact  in  Theorem  3.1,  but  first 
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we  generalize  the  definition  of  the  Bernstein  operator  to  the  range 
[a,b]. 

Definition  3.1  [15].  The  Bernstein  approximation  of  degree  n to  the 
polygon  P = [PQ , . . . » pn]  on  The  interval  [a,b]  is  given  by 

Bn[P;a,b]  = " P.(J)  (x-a)k  (b-x)n-k  (3.1) 

(b-a)n  k=0  1 k 

Theorem  3. 1 . Let  Bn [ P ; 0,1]  and  Bni[Q;l,2]  be  Bernstein  approximations 
of  degrees  n and  m to  the  polygons  P = [P^,...,Pn]  on  the  (0,1)  and 
Q = [P  , . . . ,P  ] on  (1,2),  respectively  (see  Figure  3.1).  Then 

8[8n[P]  + BJQ]]  < V[PQ Pn+in]  . (3.2) 

Proof:  (3.2)  is  easily  seen  to  be  true  if  the  sign  of  Pn  equals  the 
sign  of  P^  or  Pn+1 , since  then  V[P]  + V[Q]  = V[PQ,...  ,Pn+[T)].  Assume 
this  is  not  the  case,  i.e.,  P„  = 0 and  P , and  P are  of  opposite 
sign.  Again  (3.2)  is  true  if  either 

V[Bn[P;' 0,1 ]]  < V[P]  (3.3) 

or 

V[Bm[Q;  0,]]]  < V[Q]  , (3.4) 

since  max(V[P]  + V[Q] ) - V[PQ, . . . ,Pn+m]  + 1.  So  assume  this  is  not 
the  case,  i.e.,  V[Bn[P;  0,1]]  = V[P]  and  V[Bm[Q;  1,2]]  = V[Q] . But 
then  by  Corollary  2.9.3  Bn[P;  0,1]  and  Bm[Q ; 1,2]  must  have  the  same 
arrangement  of  sign  changes  as  P and  Q,  respectively.  Thus  (3.2) 
must  hold  in  this  case  as  well. 


I . 
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A simple  induction  argument  leads  to  the  general  statement  that 
piecewise  Bernstein  approximation  is  variation-diminishing  for  arbi- 
trarily many  segments.  A noteworthy  consequence  of  Theorem  3.1  is 
Coroll  ary  3.1.1.  Piecewise  constant  and  piecewise  linear  interpola- 
tion are  variation-diminishing. 

Indeed,  piecewise  constant  and  piecewise  linear  interpolation  corres- 
pond directly  to  B-spline  approximation  of  order  1 and  2,  respectively, 
and  the  B-spline  basis  is  known  to  be  totally  positive  of  all  orders 
[3,18]- 

Although  Bezier's  System  Unisurf  has  been  highly  successful  in 
the  design  of  automobiles,  there  are  some  problems  with  vector-valued 
piecewise  Bernstein  techniques: 

(i)  The  actual  euclidean  distance  between  the  vertices  P.. , 

i =0 , 1 m,  plays  no  role  in  the  definition  of  the  curve  segment, 

and 

(ii)  In  general,  piecewise  Bernstein  approximation  is  only  C° 
continuous.  Sometimes  design  constraints  may  require  continuity  of 
order  3 or  even  4. 

I 

Gordon  and  Riesenfeld  [19]  proposed  the  following  scheme  to  cor- 
rect for  (i):  Let  ag,  , ...  be  defined  to  be  the  fractional 

distance  of  the  i Ln  vertex  along  the  polygonal  curve  [Pg,  ...,  Pm]: 


A3 


and  define 


pi pj  ■ Bmfp5<  pr p,:t 


where  P*j  = f(j/m)  and 


f(s)  --  (ajt,  - Cj)-  [(o,jt]  - s)  Pj  + (s  - 0j)  Pjt,] 


(3.6) 


(3.7) 


for  Oj  < s < «jt). 

Note  that  (3.7)  is  just  a piecewise  linear  interpolant  to  the 


polygonal  curve  [PQ ,P-j , . . . ,P(]1]  and,  therefore,  from  Theorem  3.1  and 
Theorem  2.7  it  follows 


Cpo PJ1  i vtpo U (3-8) 

Gordon  and  Riesenfeld  observed  that,  when  the  euclidean  distances 
|Pj+1  - P j | are  all  approximately  equal,  there  will  be  little  differ- 
ence between  Bn,[po,Pl  ’ ' " and  B*[Pq,P-|  . .Pn)]*  In  extreme 

cases,  however,  the  difference  can  be  substantial,  as  shown  in 
Figure  3.2. 
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Figure  3.2  (b)  , . . . ,Pg] 

As  pointed  out  in  Section  I,  B-spline  approximation  has  been  pro- 
posed by  Riesenfeld  [5]  as  another  alternative  to  Bernstein  approxi- 
mation. Specifically,  the  B-spline  approximation  is  variation-dimin- 
ishing, has  the  convex  hull  property  and  is  a piecewise  defined  curve, 
where  the  pieces  are  joined  with  arbitrarily  high  continuity  at  the 
discretion  of  the  user.  Thus  the  B-spline  basis  has  all  the  desirable 
properties  of  the  piecewise  Bernstein  basis,  yet  without  the  inherent 
difficiency  at  the  break  points  of  the  latter  basis.  In  Theorem  3.2 
we  shall  prove  that  uniform  B-spl ine  approximation  i s variation-diminishing. 
We  first  give  a more  general  definition  of  the  uniform  B-spline  basis 
than  that  given  in  Section  I and  prove  some  elementary  properties  of 
B-spline  approximation. 

Definition  3.1.  The  ith  B-spline  basis  function  <!>.,  ...(h,x)  of  degreem-1 

~ I jIll 

with  uniform  knot  spacing  and  mesh  size  h is  defined  by 

m(h,x)  « (1/  (m-1 ) ! hm" 1 ) ) L (-l)k  ("k')  ( (mh/2)  + x - kh  - i )+1-1 
1,m  k=0  k 

(3.9) 
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where 


x > 0 

elsewhere  . 


Note  that  4* . m(h,x)  = 4.  , (h,x-h  ).  It  is  easily  seen  [8]  thatform>2 
i ) ru  1 “ 1 > ni 


. x+h/2 

4,  (h,x)  = I / <f  • ,m-l  ( h , x ) = n(h ,x)*4» . (h,x)  (3.10) 

hx-h/2  1 1 


where  * represents  convolution,  and 


7i(h,x) 


1/h  for  - h/2  <_  x <_  h/2 

0 elsewhere 


and  where 


. (h,x) 


1 for  - h/2  < x-i  < h/2 
0 elsewhere  . 


(3.11) 


(3.12) 


In  view  of  the  representation  (3.10)  we  readily  have  ^ m(h,x)  has 
positive  support  (i-hm/2,  i + hm/2)  and 


N 

Z (h,x)  = 1 , h(m-2)/2  < x < N-h(m-2)/2  (3.13) 

i=0 


where  N + 1 > nt. 

That  is,  the  convolution  of  two  positive  functions  is  positive 

and  the  finite  support  of  h = f*g,  where  f has  support  [a,b]  and  g 

has  support  [c,d],  is  given  by  [a  + c,  b + d].  As  a consequence  of 

the  local  support  of  a B-spline  basis  function,  B-spline  approximation 

is  a local  approximation  scheme.  Thus  any  finite  sum  of  the  form 
N 

F,  Pj  4>jm(h,x),  N >_  m,  involves  at  most  m nonzero  terms. 
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The  following  lemmas  will  be  integral  to  the  proof  of  Theorem  3.2. 


Lemma  3.1  [21].  Let  f(t),  { f k ( t ) , k = 1,2,., 

. . } be  real  functions 

defined  for  xe[a,b].  If 

V[fk(t)]  < n,  k = 1,2,... 

(3.14) 

and 

lim  fk(t)  = f(t), 

(3.15) 

V[f(t)]  < n. 

Proof:  Let  V[f(t)]  = N.  Then  there  exist  points  a<tQ<t^<. . .<t^b 
such  that  f(tj_-|)  and  f(tj)  are  of  opposite  sign,  j = 1,...,N.  If 
k is  sufficiently  large,  we  have  from  (3.15)  that 

sgn  fkUj)  = sgnf(tj ) j = 0,1,... ,N. 

Therefore,  for  k sufficiently  large  we  have 

n > V[fk(t)]  > V[f(t)],  (3.16) 

which  proves  the  lemma. 

Lemma  3.2.  Let  B [P:  0,n]  be  the  uniform  B-spline  approximation  of 

order  m > 2 and  mesh  size  1 to  the  polygon  P = [PQ P^]  on 

[0,n],  n >_  m.  That  is 

B"'[Pl  °'n]  = j0  Pi  *i,n.(1'x>  ' <3J7> 

where  m(l.x)  is  given  by  (3.7)  and  xe[(m-2)/2,  n-(m-2)/2]. 


47 


Then 

Zn-mZ 

VP;  O'"]  * ,p0  pj  • 

where  P™  is  defined  recursively  for  m > 2 by 

pm  B (pm-l  + pm-l)/2  ^ i = 0,1 2n_m+2 

and 


(3.18) 


(3.19) 


P? 


( P(i/2) 

^P( i-1 )/2  + P(i+l)/2)/2 


i even 

i = 0, . . . ,2n-l 

i odd 


That  is,  given  the  B-spline  approximation  of  order  m,  with  integral 
knot  spacing  to  the  polygon  P,  the  control  points  for  the  same  curve 
in  terms  of  the  B-spline  basis  over  the  refined  mesh  (0.0,  0.5,  1.0, 
. . . ,n)  are  given  by  (3.19) . 

Proof:  (by  induction  on  degree)  Let  m = 2.  We  must  show 


n 2n-m+2  ? 

E P,  p ( x)  = Z PT  ~(0. 5,x)  . (3.20) 

i=0  1=0  1 ^ 

But  the  degree  1 (order  2)  B-spline  approximation  is  the  piecewise 
linear  interpolant  to  the  vertices  P^.  Thus  (3.20)  holds,  since 
piecewise  linear  interpolation  preserves  linear  polynomials.  (See 
Figure  3.3. ) 


Figure  3.3.  Piecewise  linear  interpolation  to 

(P0 p3)  and  {Po”--’P6} 

results  in  the  same  curve. 

Now  assume  (3. 18)  holds  for  all  k,  2 £ k < n.  We  must  show  it 
also  holds  for  k = m.  From  (3.8)  we  get 

" pi  mn,x)  = l P- (-n(l  ,x)*  <fr.  ,0.x)) 

i=0  ’ i=0 

n 

= it(l,x)*  E P.  6.  ,(l,x)  (3 

1=0  ’’ 

for  Xe[(m-2)/2,n-(m-2)/2].  Now  by  our  induction  hypothesis 
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2n-ni+3 


4 Pi  V <PU-1C0-5’X)* 


thus  (3.21)  reduces  algebraically  to 


(3.22) 


2n-m+3 


E P,  4-,  _(l,x)  = tt(1,x)*  i P™'1  <)>  ,(0.5,x) 

i=0  ' i=0  ’ 


2n-m+3  , 

= (u(0.5,x-0.25)/2  +tt(0.5,x+0.25)/2)*  E P.cf>.  . (0.5,x) 

i=0  ’ 


1 2n-m+3  , 

(j)  l P-  ir(0.5,x-0.25)*d)1.jm_1(0.5,x) 


, 2n-m+3  , 

+ (£)  _E  P?'1  TT(0.5,x+0.25)*4>i>ni_1(0.5,x) 


1 / 2n -m+3  , 2n-m+3  n \ 

M i=E0  C1  *i,m(0-5’X-°-25)  + pr  ^i,m(0.5,x.0.25)j 


2n-m+3 


(3.23) 

for  xe[m-2/2,  n-(m-2)/2].  Now  both^  n^0.5,x-0.25)  and  <fr2n-m+3  m (0.5, 
x+0.25)  are  zero  on  the  interval  [(m-2)/2,  n-(m-2)/2],  and  after 

these  basic  functions  and  rearranging  the 
we  have 


i = 


dropping  terms  concerning 
remaining  terms  in  (3.23) 

Pi  fi.m'1-*’  ’ 


2n-m+2  , m , 

^ tor1  * PS*l>/2> 
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2n-m+2 

2 

i=0 


P,n 

i l ,m 


(0. 5 ,x) 


(3.24) 


Thus,  we  have  proved  the  lemma. 

n 

Lemma  3.3.  Let  B [P;  0,n]  = E P.  <p.  (1 , x ) as  in  (3.17)  and 

^ _Q  I l jlTI 

define 


*0[P]  - CPS.»T 


where  P™  is  given  by  (3.18)  and 


*k[P]  - *0[/"’[P]]  > [P^'k P”('kJ  , (3.26) 

for  p(m)  = 2(2n-m+2)k  - m+1 , k ^ 1 . 

Then, 


lim  K»k[P]  = Bm[P;  0,n]  . (3.27) 

k-HO 

Riesenfeld  [33]  has  recently  given  a proof  of  Lemma  3.3  for  the 
case  m=3.  Here  we  use  a quite  different  approach  for  the  proof  for 
arbitrary  but  finite  m. 

Proof:  By  Lemma  3.2  and  induction  on  k we  have 

Bm[^k[P];  0,n]  = Bm[P;  0,n]  . 


Mow  by  the  minimizing  the  nature  of  piecewise  linear  interpolation  it 
is  easily  seen  that 
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arc  length  (/[P])  < arc  length  (^k_1[P])  . (3.28) 

By  applying  the  triangle  inequality  to  (3.18)  we  can  deduce 


iCi  - 0 i’/2  <arc  len9th  cC' pp(',-i)i> 

which,  combined  with  (3.28)  and  induction,  imply 

|pm,k  _ pm,kj  £ (1/2m+k)  L>  (3.29) 

where  L is  the  arc  length  of  [P0,...,Pn]  and  | | means  euclidean 
distance.  That  is,  by  (3.29) 


But  then 


lim  |P*;k  - Pn.',k|  = o 

k-«o 


(3.30) 


lim  |pT;k  - Pm’k 
k-K-  1+J  1 


(3.31) 


for  arbitrary  ieO  , . . . ,p(m) } , j = l,...,m. 

From  (3.27)  and  (3.13)  we  know  Bm[P ; 0,n](xg),  for  each 

xQc( (m-2)/2,n-(m-2)/2) , lies  within  the  convex  hull  of  P'!1,k,  PI!1jk, 

ro  k 

".’Pi+in  ^or  some  w1*""  (3-31),  we  can  then  conclude 

lim  *k[P]  = Bm[P;  0,n]. 

k-HO 


(3.32) 


52 


Theorem  3.2.  Uniform  B-spline  approximation  is  variation-diminishing. 
Proof:  We  have  from  (3.32) 

lim^k[P]  = B [P;  0,n]. 

k-«o  m 

It  follows  from  the  piecewise  linear  nature  of  the  construction  (3.19) 
and  Corollary  3.1.1  that  V[^[P]]  £ V [ P ] for  all  k.  The  conclusion 
now  follows  from  Lemma  3.1. 

Bernstein-Bezier , B-spline  and  other  generalizations  of  the  Bezier 
curve  [13,14]  still  appear  to  have  inherent  shortcomings  for  realtime 
interactive  design.  In  particular, 

I.  It  is  often  the  case  that  a user  wishes  to  create  a local 
fit  to  the  polygon  in  his  design,  yet  there  are  no  'natural' 
handles  or  control  parameters  in  the  above  methods  for  the 
designer  to  manipulate  interactively  in  order  to  achieve 
these  shapes. 

II.  Further,  the  extension  from  curves  to  surfaces  with  the  above 
schemes  requires  a rectilinear  network  of  control  points,  a 
severe  restriction  on  interactive  design. 

For  instance,  the  addition  (deletion  of  a point  to  a rectilinear 
net  requires,  in  general,  the  addition  (deletion)  of  a whole  row 
and  column  of  points  in  order  to  retain  the  rectilinear  topology 
(see  Figure  3.4). 
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V(7,7) 


Figure  3.4.  Rectilinear  topology.  Note  that  to  remove 
the  center  vertex  and  still  retain  the 
rectilinear  structure  one  has  to  remove  all 
of  Row  4 and  Column  4. 


The  Construction  of  New  Models 


We  now  apply  the  theory  of  total  positivity  to  the  construction 
of  new  models  for  curve  and  surface  design,  with  emphasis  on  methods 
which  avoid  the  deficiencies  I and  II  discussed  in  the  previous  sec- 
tion. Specifically,  for  curves  we  develop  linear  operators 


L[P;a,b]  = l P^-U.c^)  for  P = [PQ, . . . .Pj  tc(a,b) 


(3.33) 
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where  the  can  be  varied  by  the  user  of  the  system  to  control 
local  "closeness"  of  fit  to  the  polygonal  curve  P,  thus  giving  the 
user  the  ability  to  create  "cusps"  and  "flats"  in  his  design  with 
the  same  natural  flexibility  he  has  in  moving  the  vertices  of  the 
polygon.  Further,  we  build  in  the  desirable  properties  of  the 
Bernstein  and  B-spline  methods, 

V(L)  <_  V[P]  (variation-diminishing  property)  (3.34) 

and 

(t,a. ) >_  0, 

for  all  i with 

E <{>.j(t,a.j)  = 1 (convex  hull  property)  (3.35) 

It  is  well  known  [15,16]  that  B-spline  approximation  to  a continu- 
our  function  f on  an  interval  [a,b]  for  fixed  degree  converges  to  f 
as  the  mesh  size  h goes  to  zero.  In  terms  of  our  primitive  polygon 
P,  we  can  get  a closer  local  fit  to  the  polygonal  curve  P if  we 
sample  not  only  at  the  vertices  of  the  polygon  but  also  in  the  inter- 
val in  which  we  wish  to  approximate  more  closely.  The  more  samples 
in  the  region  of  interest  the  closer  our  approximation  to  the  polygon 
there.  That  is,  if  we  let  the  polygon  P be  defined  by 

P(t)  = E P.  <j>.(t)  (3.36) 

j=0  3 3 

where  the  <j> . are  the  piecewise  linear  cardinal  functions.  Rather 

J 

than  form  the  B-spline  approximation  to  the  polygon  P by  sampling  at 
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the  vertices  P.,  we  define 

J 


G(t)  " kf,  Wt)p(tk> 


(3.37) 


where  the  4>,  are  the  B-spline  bas-is  functions  of  degree  m and  the 

K * ill 

1 ^ are  the  knots,  which  are  not  just  located  at  the  vertices,  but  in 
clusters  in  intervals  of  local  interest  (see  Figure  3.5). 


Since 


p(tk>  ” ^ pj  ♦A*  • 


(3.37)  can  be  rewritten  as 


G(t)  = E <C(t)  P(tJ 
k=l  k K 


l *!"(t)  ( l P,  <Mt.  )) 
k=l  K j=0  J 3 K 


n p 


E PA  l C( t)  *.(t,)) 
j=0  J k=0  K J K 


E P*  +•  m(t)  , 


J J 


where 


'j.JO  ■ * *k(0  +jCtk) 
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Since  the  B-spline  basis  is  TP,  we  know  from  Theorem  2.6  and 

Theorem  2.7  that  the  basis  {<t>i  } is  TP.  Thus  by  appropriately 

choosing  the  knots  { t.  } we  can  generate  a new  basis  {<{>'.}  which  more 

K J 

closely  approximates  the  polygon  P,  yet  retains  the  valuable  convex 

hull  and  variation-diminishing  properties.  More  importantly,  this 

construction  method  generalizes  to  other  TP  bases: 
n 

0 ( x ) is  TPn  on 

XxJ,  J = {l,2,...,n}  and  define  recursively 


Let  Gq(x)  = E if.  0(x)  P.,  xe[a ,b] , where  <#> . 

n=i  J>  J J 


G^x) 


n ( i ) 

^ ^ (x)  G.  I (x  • • ) , 

j=l  1-1  J»' 


i = 1 ,2 , . . . ,m  (3.38) 


where  xe[a,b]  and  <t»j  i (x)  is  TP^.j  on  XxJ(i),  J(i)  = {1  ,2 , . . . ,n(  i ) } 
where  n(i)  _>  n(i-l)  and  the  sequence  (x-  .}  is  the  ith  knot  vector. 

J J ' 

Expanding  G.  ,(x.  • ),  (3.38)  can  be  rewritten 

* “ ' J > ' 


Gi(x) 


n(i) 


n(i-l ) 


Z i(x>  1 i-l(xi  i1  Gi  2(xk  i 
j=l  J>1  k=l  k)1  1 1 1 k>1 


n(  i-1 ) / n ( i ) 


n(i-l) 

kf1  Gi-2^xk,i-l)  ^k,i-l(x^ 


(3.39) 


By  Theorem  2.7  <f.£  ^(x)  is  anc*’  therefore,  by  Theorem  2.10 
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< V[Gj.2(xjij.,),..,Gi.2(xn(i.2)>i.1)]  (3.40) 


Repeating  the  above  algorithm  we  readily  deduce 


G (X)  = r.  d>'.  (x)  p., 

1 j=l  J,u 


(3.41) 


where 


V[G.(x)]  £ V[Pr  Pn] 


(3.42) 


Now  let 

n 

<}>(x)  = !<}>•  n(x) 
j=l 

Then 

<t>':  n(x)  = 4>-  _(x )/<J>(x)  , j = l,...,n 

J J >u 


is  TP  by  Theorem  2.4  and 

n J 


£ 4>'l  0(x)  = 1 

3=1 


(3.43) 


(3.44) 


It  will  now  be  demonstrated  how  the  techniques  developed  in  (3.33) 
and  (3.43)  for  the  construction  of  the  variation-diminishing  and 
convex  hull  properties,  respectively,  can  be  ust^  to  generate  a set 
of  basis  functions  with  the  "natural"  handles  desired.  Let 


G-|  (x)  = P = [P1 Pn]  (3.45) 

be  the  piecewise  linear  interpolant  to  the  points  { P ^ , i = l,...,n) 
over  the  uniform  mesh  {l,...,n}.  Define 


60 


Now  setting  c»j  1 - «j_-j  • for  j odd,  and  rearranging  terms  as  in 
(3.38)  we  have 

n 

G2(x)  = .j,  .j(|)  Pj,  (3.48) 

where  the  aj  ^ (x'»aj_i  >aj  ,aj+i ) are  TP  ar|d  vary  in  shape  from  trape- 
zoidal (ctj_i  ,ctj  >ctj+i  = 0)  to  triangular  (aj_1  ,a.  ,a.+}  = 1)  (see 
Figure  3.7). 

Now  let  {ipj(x),  j = i , . . . ,n } be  any  set  of  totally  positive  functions 
and  define 


(3.49) 

(3.50) 

For  example,  let  ipj (x)  be  the  uniform  B-spline  basis  of  degree  m-1 . 
Then  for  a.  small,  the  basis  functions  are  "bell -shaped, " while  for 
large,  these  basis  functions  converge  in  shape  to  the  piecewise 
linear  basis  functions.  This  convergence  takes  markedly  for  the 
lower  degrees  (see  Figures  3.8  and  3.9). 


G (x)  = E iK(x)  G?(x  • ) 

o j=l  J * J 


' * *j<x>  f Vr“i"W  h 

J ' J • 


i = l 


^x;ai-l  >ai»ai+-| ) Pi > 


where 


V-j  (x,u-j_i  >aj » + 1 ) ^ 'l'j(x)  (J),-  t (x^ ; otj  i ,a- ) 

0 • 


T ' >1 ' j i-1  i i+1 


Figure  3.9  (a).  Canonical  basis  function  4'q ( x ; 0 , 0 , 0 ) for 
the  choice  ip  (x)  = B-spline  basis  of 
degree  5. 
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The  previous  construction  give  us  an  indication  of  the  possibili- 
ties for  creating  new  TP  bases  from  other  TP  bases.  If  in  (3.49)  we 
choose 

n-l  . , . 

i|k(x)  = £ (n.  ) x^l-x)"  1 , i.e.,  the  Bernstein  basis  of  degree 

J i=0  1 

n-l,  then  our  construction  is  very  similar  to  that  of  Gordon  and 
Riesenfeld  defined  by  equations  (3.6)  and  (3.7),  only  here  we  have 
allowed  the  number  and  position  o the  Pj  to  be  flexible. 

There  are  other  schemes  which  easily  can  be  developed  from  a 
careful  perusal  of  the  construction  prototype  (3.38).  For  instance, 
if  in  (3.49)  we  take  {ip  - (x) , j=l , . . . ,n}  to  be  the  B-spline  basis  of 

\J 

degree  2 and 

ci j — 1/2,  j * 1 , . . . ,n , 

then  the  construction  reduces  to  that  of  (3.18).  From  Lemma  2.5 
we  know  that 

G3(x)  = B3[P;  0,n]. 

The  method  can  be  extended  by  generalizing  the  construction  (3.19) 
to  create  a whole  family  of  bases  n)}  such  that 

(x;l/2,l/2,l/2)  = <fj>m(x),  j = 1 n,  (3.51) 

where  . (x)  is  the  mth  order  B-sp.1  i ne  basis  functions, 

j ,in 
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Thus  the  family  of  bases  {<!<'.  , i = l,...,n}  is  a proper  gener- 

I >111 

alization  of  the  family  of  B-spline  bases,  where  the  {>j/  m>  i = 1,..., 
n}  form  a substantial  improvement  in  terms  of  the  'natural'  handles 
desired.  (See  Figures  3.10  and  3.11.1 

More  Curve  Techniques 

In  the  previous  constructions  we  used  a discrete  form  of  Theorem 
2.7  to  construct  TP  functions  from  other  TP  functions.  An  equally 
valuable  technique  for  CAGD  results  when  the  variable  of  summation 
in  (2.16)  is,  in  fact,  continuous.  For  example,  let 

n 

P(t)  = E <t>,(x)  P.  (3.52) 

i=0 

as  defined  in  (1.1),  and  define 

CO 

Q(t)  = T(P) ( t)  = / K(t,s)  P(s)  ds,  (3.53) 

-00 

2 

where  K(t,s)  is  TP  on  R . Then  by  Theorem  2.2  we  have 

V[Q]  < V[P]  . (3.54) 

Rewriting  (3.53)  we  have 

Q(t)  = / K(t,s)  P(s)  ds 

-00 

= / K(t,s)  E ^(s)  Pi  ds 

-oo  i=0 
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n °°  ' ' 

= E (/  K(t,s)cJ).(s))  ds  P. 

i=0  1 

n 00 

= e 4.1  (t)  P.,  where  4>l(t)  = / K(t,s)<}>.  (s)  ds 

1=0  ' -CO 

In  view  of  Theorem  2.7,  provided  that  [4>^  ( t ) 3 is  TPn>  it  follows 
C0j(t)]  is  TPn. 

To  develop  a class  of  kernels  K(t,s)  suitable  for  CAGD  we  re- 
strict our  attention  to  TP  kernels  of  the  form  K(x,y)  = f(x-y) 

2 

(translation  kernels)  where  (x,y)  e R and  / f(x)  = 1.  Under  these 

-oo 

constraints  the  transformation  (3.53)  becomes  a convolution,  i.e., 


T(P)(t)  = / K(t-s)  P(s)  ds.  (3.56) 

-CO 

Recall  from  Section  II  that  the  function  exy  is  TP  on  R*\  Then 
by  Theorems  2.4  and  2.5, 

K,(x,y)  = e-“(«-y)2 

_ e-ax2  g-ay2  g2axy  (3.57) 

is  TP  for  a > 0.  Similarly,  it  can  be  shown  that 

K2(x,y)  = e'a|x_yl  (3,58) 

2 

is  TP,  for  a > 0,  (x,y)  e R . Since  K^(x,y)  and  K2(x,y)  are  TP, 
from  Theorem  2.7  we  know  their  self-convolutes  are  TP.  That  is, 
if  we  define 
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i 

I 

K^(x,y)  = K^fx.y)  * K.(x,y)  i = 1,2;  j = 1,2..., 

(3.59) 

where 

l 

K°(x,y)  = K.(x,y)  i = 1,2, 

then 

K^(x,y)  is  TP  on  R2.  (3.60) 

The  kernels  K'j  are  "bel  1 -shaped"  and  symmetric,  the  "spread"  of  the 
curves  depending  on  a.  (See  Figures  3.12,  3.13.)  Note  that  K'j(x.y) 
is  of  continuity  class  c”  for  all  j,  while  K^x.y)  e CJ. 

If  P(t)  in  (3.52)  is  our  primitive  polygon,  which  means  that 
(<Mx)}  is  the  cardinal  piecewise  linear  basis,  then  we  can  construct 
new  bases 

4> j ( x)  = Kjj(x)  * 4>i  (x)  (3.61) 

for  p = 1,2;  i = 0,1,... ,n;  and  j = 1,2 such  that  {^(x)}  is 

TP  on  Xxl,  X = (-“,«>),  I = (0,1,.  ..,n}  for  all  j. 

There  are  several  alternative  methods  of  ensuring  that  the  bases 

j 

{^(x))  enjoy  the  convex  hull  property.  We  can  either  use  the  tech- 
nique developed  in  (3.43)  or  we  can  normalize  (K'j(x)}  such  that 

00 

/ k](x)  dx  = 1. 

.00 
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That  is,  define 

L'j(x)  = K-j(x)/C  , i = 1,2;  j = 1,2,  ... 

where 

oo  m 

C = / K;j(x)  dx. 

-OO 

Then  from  (3.61)  we  have  for  i = 1,2, 

n i n . 

E 4(x)  = E (q(x)  * <f»k(x) ) 
k=0  k=0  1 

i n 

= M(x)  * E 4>k(x) 

1 k=0 

(3. 

= L^(x)  * 1 

= 1 

Of  course,  in  view  of  the  recursive  nature  of  (K'j(x,y)},  various 
combinations  of  (3.43)  and  (3.62)  could  be  used  as  well. 

Net  surprisingly,  in  view  of  their  similarity  in  shape  to  the 
B-spline  basis  (compare  Figures  3.14  and  3.15  and  Figures  3.16  and 
3.17),  curves  formed  with  these  new  bases  can  be  remarkably  close  in 
shape  to  the  B-spline  approximation.  This  similarity  is  explained 
mathematically  by  the  Central  Limit  Theorem. 

Theorem  3.3  (Central  Limit  Theorem)  [27].  Let  f^(t),  t (-°°,»),  i 
an  integer,  be  real  positive  and  symmetric  about  t = 0,  such  that 

/ ft(t)  dt  = 1 
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for  all  i and  define 

f(t)  = f,(t)  * f2(t)  fn(t),  (3.63) 

where,  again,  * represents  convolution.  Then  if 
“ ■» 

/ tJ  f.(t)  dt  < C (3.64) 

CO 

where  C is  an  arbitrary  constant,  and 

n «> 

lim  E / rfi(t)dt=“  (3.65) 

n-*°  i = l -°° 

Then 

lim  f(t)  = — — e-(t)2/2o2  (3.66) 

n-xo  o/2tT 

where 

a2  = E / t2f.(t)  dt 
i=l  -»  1 

We  constructed  the  {rj^(x)}  in  (3.61)  precisely  as  the  iterated 
convolution  of  j functions,  and  a reexamination  of  the  B-spline  basis 
of  degree  m reveals  that  each  is  the  iterated  convolution  of  m func- 
tions, as  well.  Without  going  into  detail,  it  can  be  shown  [3,12] 
that  all  of  these  functions  satisfy  the  hypothesis  of  Theorem  3.3. 
Therefore,  in  the  limit,  all  these  bases  would  be  of  Gaussian  form, 
differing  only  in  their  dispersion  about  their  point  of  symmetry. 

Thus  the  similarity  in  shape  of  the  various  approximation  methods 
is  in  actuality  a reflection  of  this  "tendency  to  Gaussian  form"  in 
the  basis  functions. 


75 


Surfaces 


Although  the  extension  of  vector-valued  curve  methods  to  surfaces 


is  straightforward,  there  is  no  satisfactory  theory  of  total  positiv- 
ity for  functions  of  more  than  one  variable  [3].  However,  when 
dealing  with  surface  equations  of  the  form 


n 

S <Mx,y)  P.  (3.67) 

i=0 


we  can  refer  to  the  total  positivity  of  the  {<t> • (x ,y ) > with  respect  to 
x and  i or  y and  i,  respectively.  As  in  our  development  for  curves, 
we  are  interested  in  bases  {<j»i(x,y)}  which  are  totally  positive  with 
respect  to  both  continuous  variables. 


Let  L-|  and  be  some  linear  operators  over  the  polygonal  curves 
P = [P0i-...Pn]  and  Q = [Q  ,...,Qn],  respectively,  defined  by 


n 

L,[P]  = I d>- (x)  P.,  xeR 

i=0  1 

and 

m 

L?[Q]  = Z ip.(y)  Q. , yeR 
c i=0  1 1 


where  the  bases  0 (x),  i = 0,...,n]  and  [ij'.(y),  j = 0,...,m] 
i J 

are  TP  on  their  respective  domains  and  possess  the  convex  hull  pro- 
perty. Then  given  a rectilinear  network  R = [P. ;,  i = 0,...,n, 
j = 0,...,m]  we  can  define 


n 

L,[RJ  = Z 4> - ( x ) P..  for  j fixed,  jc  [0,  m]  (3.68) 
i=0  J 
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L?[R]  = E <J/.(y)  P..,  for  i fixed,  ic  [0,  ....  n]  . (3.69) 

j=0  J J 

L1  and  are  called  lofting  operators,  reflecting  the  fact  that  they 

filter  or  vary  the  shape  of  R with  respect  to  only  one  variable.  In 
view  of  the  total  positivity  of  [<fc.(x)]  and  0-(y)],  we  have 

V(L1[R]](x)  < V [ R ] ( i ) 
and 

V(L2[R]](y)  < V[R]( j) , 

where  we  have  altered  our  original  notation  for  the  number  of  sign 
changes  for  curves  to  reflect  which  domain  we  are  considering.  If  we 
wish  to  smooth  in  both  directions,  we  construct  the  tensor  product 
of  the  two  lofting  operators  as 

L[R]  - L1 [L2[R]]  = L2[L1[R]],  (3.70) 

where  we  have 

V[LrR]](x)  < V[L1 [R]](i ) 

V[L[R]](y)  < V[L2[R]](j) 

Although  these  kinds  of  generalizations  from  curves  to  surfaces 
have  proven  extremely  successful  for  CAGD  [121,  and  we  are  dealing 
formally  with  two-dimensional  surfaces,  the  approach  is  inherently 
one-dimensional.  This  fact  is  reflected  in  the  polyhedral  network 
for  the  lofting  operators,  and  therefore  in  the  tensor  product  oper- 
ators, which  manifestly  require  rectilinear  control  points. 
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There  is  a curve  to  surface  generalization,  modelled  after  the 
convolutional  methods  for  curves  (3.56),  which  circumvents  the 
use  of  rectilinear  networks.  Let 


P = l <|>,(x,y)P.  (x,y)cR  (3.71) 

1=0  1 

3 

where  the  P.  are  arbitrary  points  of  R corresponding  to  knots 
2 

= (x^y^)  of  R and  where  the  {<}>.}  are  piecewise  linear  poly- 
nomials of  two  variables  which  solve  the  interpolation  problem. 


( 1 i = j 
l 0 elsewhere, 


0 < i , j <.  n . 


Thus,  for  distinct  P^  P is  a proper  triangul ation  of  the  points  P^ 
(See  Figure  3.18.)  It  is  evident  that  a given  polygonally  bounded 
domain  in  the  plane  can  have  several  triangulations  (Figure  3.19). 
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For  various  efficient  algorithms  for  triangulation  of  the  plane,  see 
[23,24]. 

o 

Now  for  K(x,y)  and  M(x,y)  TP  on  R as  in  (3.56),  where  K(x,y) 
and  M(x,y)  are  translation  kernels,  define  the  lofting  operators 


and 


L-j  [P] 


CO 


I 


n 

K(x,s)  ( Z <J> - ( s,y ) P,)  ds 
i=0 


L2[P] 


/ 


n 

M ( x , t ) ( Z <J>.  (x,t)  P.)  dt 
i=0  1 1 


As  in  (3.68)  and  (3.69)  we  have 


(3.72) 


(3.73) 


VtL^PJKx)  < V [P ] ( x ) 
and 

V[L2[P]](y)  < V[P](y)  . 

In  view  of  the  piecewise  linear  nature  of  the  we 

have 

oo 

<t>-(x,y)  = }m  K(x,s)  4>-(s,y)  ds  (3.74) 

is  TP  on  Xxl , and 

on 

4>-(x,y)  = / M(y,t)  <f>.(x,t)  dt  (3.75) 

-oo 

is  TP  on  Yxl,  where  XxY  = and  I = (0,1, ...,n}. 

Of  course,  we  can  generalize  the  lofting  transformations  to  the 
tensor  product  operation  by 
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L[P]  = L-,[L2[P]]  = L2  [L  i [P]  ] 

Although  we  have  insisted  in  (3.70)  that  the  <j>.  represent  the 
cardinal  piecewise  linear  basis  functions,  the  development  which 
follows  (3.70)  works  equally  well  if  the  represents  the  piecewise 
bilinear  basis  functions  associated  with  the  vertices  of  a rectil inear 
network.  Thus  our  technique  of  approximating  triangular  networks 
of  points  is  a proper  generalization  of  tensor  product  approximation 
to  rectilinear  networks,  encompassing  the  latter  as  a special  case. 
Note  that  in  avoiding  the  biased  directions  of  approximation  inherent 
in  the  rectilinear  schemes,  we  must  choose  the  directions  of  approxi- 
mation. Above  we  chose  to  integrate  with  respect  to  the  x and  y 
directions,  but  clearly  any  two  independent  directions  would  have 
sufficed.  The  following  figures  show  various  tensor  product  approxi- 
mations to  their  respective  triangular  nets  for  different  choices 
for  L and  M. 
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Figure  3.22. 


T * (e 


-We-|y|,2. 


Figure  3.23.  T * (e'lxU‘lyl)4. 


IV.  CONCLUSION 
Summary 

The  techniques  developed  and  extended  herein  for  curve  and  sur- 
face design  form  a proper  generalization  of  Bernstein-Bezier  and  B- 
spline  methods,  while  enjoying  increased  flexibility  for  interactive 
manipulation. 

The  application  of  the  theory  of  total  positivity  to  CAGD  has 
not  only  resulted  in  well  structured  approaches  to  new  mathematical 
modelling  for  ab  initio  design,  but  a new  framework  for  analyzing 
and  understanding  existing  techniques. 

Although  this  paper  has  not  directly  attacked  the  problems  of 
computability  for  the  new  bases,  the  stability  of  the  methods  is 
inherent  in  the  "bell-shaped"  form  of  the  basis  functions,  as  exem- 
plified in  their  tendency  to  Gaussian  form.  Further,  the  construction 
(3.18)  forms  an  efficient  algorithm  for  calculating  B-splines  and  can 
be  extended  to  the  other  bases  constructed  in  Section  m.  These 
construction  methods  form  the  core  for  a class  of  geometric  algorithms 
for  computing  the  derivative,  arc  length  and  intersections  of  spline 
curves,  as  well  as  the  area,  volume  and  intersections  of  the  corres- 
ponding spline  surfaces. 

The  feasibility  and  utility  of  the  newly  constructed  models  for 
interactive  design  are  demonstrated  in  the  following  figures,  which 
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are  "frames"  from  an  interactive  session  on  an  experimental  system 
for  curve  and  surface  design.  Since  the  final  decision  on  what 
constitutes  a "good"  design  is  subjective  and  will  probably  vary 
from  user  to  user,  it  is  felt  that  the  ability  developed  herein  to 
vary  the  mathematical  model  while  retaining  the  desirable  variation- 
diminishing  and  convex  hull  properties  may  be  of  increased  impor- 
tance in  future  CAGD  systems. 

Overview 

The  authors  feel  that  this  constitutes  a natural  third  paper 
in  the  direction  of  establishing  a sound  mathematical  basis  for  the 
application  of  the  variation  diminishing  approximation  method  to 
computer  aided  geometric  design,  the  first  two  papers  being  Gordon 
and  Riesenfeld  [5,28]. 


Figure  4.3  (a).  Fifth  degree  B-spline  approximation  to  the 


polygon  given. 


Figure  4.3  (b).  Fifth  degree  spline  approximation  with  the 


basis  constructed  from  (3.50). 
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Figure  4.4  (e).  A triangulation  T of  S. 
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Figure  4.5  (b).  Cubic  B-spline  approximation  to  S.  Note 
there  are  only  local  differences  in  shape 
from  Figure  4.4  (b). 


Figure  4.5  (c).  A triangulation  T of  S. 
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2 2 

Figure  4.5  (d),  T * (3/2tt)  e"3x  ^ e"3y  /'2. 

Note  there  are  only  local  differ- 
ences in  shape  from  Figure  4.4  (f). 
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